function [err] = depfits2(varargin)
% depression/recovery fits of synapse data.
% This version is for single cells and multiple frequencies...
% It operates directly off the CONTROL table after analysis has been done
% Call with no inputs; uses the selected data sets to build the tables we
% need in order to do the fits. Frequency is extracted from the parameter
% field... just as in the analysis routine (see EPSC_train.m).
% If this is called with an argument list, the first argument is the
% paremeter set for the initial fits.
% 9/26/2008 Paul B. Manis
%

global CONTROL

sf = getmainselection();

optimize = 1; % set this flag to do the fitting; otherwise we just plot the data.
plota = 0;

clist = {'ro-', 'ys-', 'g^-', 'b+-', 'c.-'};

if(nargin == 0)
    % synapse parameters:
    % (Dittman and Regehr model)
    table.F = 0.4; % release probability (constant; no facilitiation)
    table.k0 = 1/1.75; % /s, baseline recovery rate from depletion (slow rate)
    table.kmax = 1/0.025; % /s, maximal recovery rate from depletion (fast rate)

    table.td = 0.05; %  time constant for calcium-dependent recovery
    table.kd =  0.7; % affinity of fast recovery proces for calcium sensor

    table.ts = 0.015; % decay time constant of glutatme clearance
    table.ks = 1000; % affinity of receptor desensitization for glutatmate
    % The large value means no desense occurs (otherwise, ks should be about
    % 0.6)

    table.kf = 0.6; % affinity of facilitiation process
    table.tf = 0.01; % make facilitation VERY slow

    table.dD = 0.02; % sets Ca that drives recovery(Ca influx per AP)
    % 0.02 yields rate-dep recovery in 100-300 Hz
    table.dF = 0.02; % sets Ca that drives facilitation

    table.glu = 0.3;
else
    p = varargin{1};
    % note: order is important!!!!!
    table.F = p(1);
    table.k0 = p(2);
    table.kmax = p(3);
    table.td = p(4);
    table.kd = p(5);
    table.ts = p(6);
    table.ks = p(7);
    table.kf = p(8);
    table.tf = p(9);
    table.dD = p(10);
    table.dF = p(11);
    table.glu = p(12);

end;

j = 1;
freq = zeros(length(sf), 1);
for i = sf % get the data sets we need.
    par =  make_struct(CONTROL(i).Parameters, 'delay 5, np 20, freq 20, tau_init 40, nrep 4');
    freq(j) = CONTROL(i).EPSC_Train.Freq;
    tb{j} = CONTROL(i).EPSC_Train.stim;
    ib{j} = mean(CONTROL(i).EPSC_Train.amps, 1);
    trb{j} = CONTROL(i).EPSC_Recovery.Tdelays-CONTROL(i).EPSC_Recovery.rdelay;
    rb{j} = mean(CONTROL(i).EPSC_Recovery.Iamps, 2)';
    j = j + 1;
end;
%trb={};
%rb={};
[fs, forder] = sort(freq);



for i = 1:length(tb)
    tb{i} = 0.001*tb{i}; % conversion to seconds
end;
for i = 1:length(trb)
    trb{i} = 0.001*trb{i}; % conversion to seconds
end;

newfigure('depfits', 'Depr. Fits');
if(plota)
    subplot(3,1,1);
    for i = 1:length(tb)
        plot(tb{i}, ib{i}, clist{forder(i)});
        hold on;
    end;
    subplot(3,1,2)
    for i = 1:length(trb)
        hr = plot(trb{i}, rb{i}, clist{forder(i)});
        hold on;
    end;
    set(gca, 'XScale', 'log');
end;

err = 0;
if(plota)
    subplot(3,1,3)
end;
for i = 1:length(tb)
    if(~isempty(trb))
        ntb{i} = [tb{i} max(tb{i})+trb{i}];
        nr{i} = [ib{i} rb{i}];
    else
        ntb{i} = [tb{i} tb{i}(end)];
        nr{i} = [ib{i} ib{i}(end)];
    end;
    nr{i} = nr{i}/nr{i}(1);
    plot(ntb{i}, nr{i}, clist{forder(i)}(1:2), 'MarkerFaceColor', clist{forder(i)}(1), ...
        'MarkerSize', 2.25);
    hold on;
    %    set(gca, 'Xlim', [0 0.5]);
    %    set(gca, 'Ylim', [0 2.0]);

end;

if(optimize)
    [xo, yo] = XuF(2, table, tb, trb);
    for i = 1:length(xo)
        hold on
        plot(xo{i}, yo{i}, clist{forder(i)}, 'markersize', 0.1);
    end;

    err = 0;
    for i = 1:length(xo)
        err = err + sum((yo{i}' - nr{i}).^2);
    end;

end;





















